Quantum Monte Carlo study of the three-dimensional attractive Hubbard model 
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We study the three-dimensional (3D) attractive Hubbard model by means of the Determinant 
Quantum Monte Carlo method. This model is a prototype for the description of the smooth crossover 
between BCS superconductivity and Bose-Einstein condensation. By detailed finite-size scaling 
we extract the finite-temperature phase diagram of the model. In particular, we interpret the 
observed behavior according to a scenario of two fundamental temperature scales; T* associated 
with Cooper pair formation and with condensation (giving rise to long-range superconducting 
order). Our results also indicate the presence of a recently conjectured phase transition hidden by 
the superconducting state. A comparison with the 2D case is briefiy discussed, given its relevance 
for the physics of high- Tc cuprate superconductors. 



The existence of a smooth crossover between the 
two paradigms of quantum superfluidity, the Bardeen- 
Cooper-SchriefFer (BCS) superconductivity and the 
Bose-Einstein condensation (BEC) is firmly established 
[0, ||. In this context, the attractive Hubbard model 
(AHM) has appeared as an ideal presentation of the 
whole evolution between the BCS and BEC physics |^ . A 
concrete property of this Hamiltonian is the existence of 
two (not always) distinct energy scales: one associated 
with the formation of Cooper pairs (T*) and another 
with the onset of long-range order in the system (Tc) 
[Q. Although their qualitative behavior is well-known, 
a quantitative determination is still missing, due to the 
fact that it is hard to access the intermediate regime by a 
controlled approximation scheme. In this respect the De- 
terminant Quantum Monte Carlo (DQMC) method |^ 
is a powerful tool as it provides results free of systematic 
errors. A detailed finite-size analysis is however necessary 
in order to extract the thermodynamic limit properties, 
which can then be compared with the outputs of other 
methods recently applied to the same problem @, §. At 
this point we should stress the role of dimensionality that 
determines the nature of the superconducting phase tran- 
sition at Tc, the strictly 2D realization of the model is 
characterized by a Berezinskii-Kosterlitz-Thouless type 
phase transition, whereas the 3D case displays a "nor- 
mal" second-order one, which is more easily accessible by 
DQMC. Since the intermediate regime of the AHM con- 
stitutes the simplest model for a short-coherence-length 
superconductor, the considerations presented hereafter 
may as well help to clarify the influence of the dimen- 
sionality on some properties exhibited by the 3D strongly 
anisotropic high- Tc superconductors. 
In this Letter, we present the results of extensive DQMC 
simulations for the finite-temperature properties of the 
AHM in three dimensions. In spite of finite-size ef- 
fects, we show that it is possible by a scaling analysis 
to quantitatively establish the phase diagram of Tc(C/, n) 
as a function of the interaction strength and density of a 
model that exhibits a genuine second-order phase tran- 



sition (unlike its 2D version). Furthermore, the pair for- 
mation temperature T* is studied in detail, revealing the 
existence of a transition in the non-superconducting state 
taking place at a critical coupling strength. These results 
complete recent calculations which have postulated the 
existence of such a transition in the infinite-dimension 
version of the model || . 

Model and method. — The attractive Hubbard model is 
defined by the following Hamiltonian, 
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where < i,j > denotes a pair of nearest neighbors on a 
cubic lattice with N = sites, c\^{cicr) is a fermion cre- 
ation (annihilation) operator of spin a =|| and Ui^a = 
cl^Cia- We take t > 0, U > and the chemical potential fi 
is tuned to yield a fixed density < n < 2. Outside half- 
filling (n 7^ 1) this model presents a finite-temperature 
transition into a phase characterized by long-range s- 
wave superconducting order associated with the breaking 
of the U(l) gauge symmetry. 

To study the finite temperature properties of this sys- 
tem we use the conventional DQMC H, |^ simulation 
method. Since the attractive interaction does not lead to 
a minus-sign problem, the whole U-n-T phase diagram 
can be reliably studied. Because of the grand-canonical 
nature of DQMC, it is necessary to estimate the function 
fi — fJ.{T, n, U, L) in order to work at a fixed density n. 
This presents a considerable load in this work compared 
to similar DQMC simulations at half-filling |^ . Typically 
we take n = 0.5 (quarter-filling) for which results using 
other methods have already been presented ||^, ||]. We 
also restrict ourselves to finite-temperature static corre- 
lation functions such as the s-wave pair-pair corre- 
lation function Ca and the Pauli spin susceptibility xp. 
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Here A^ = Ci|Cix and S,; = I]^,i/=t,i '^L'^'^'''^^''' being 
the vector of Pauli matrices. Ca allows to determine 
the superconductiong transition temperature Tc, since it 
signals the breaking of the U(l) gauge symmetry. We 
recall that this approach is not applicable to the strictly 
2D case where more sophisticated quantities have to be 
calculated On the other hand xp indicates the pres- 
ence of pairing in the system, related to the temperature 
scale T* as discussed below. 

Regarding the DQMC simulations, the imaginary time 
discretization is At = 0.125i^^ and lattices of size 
N = — IQ^ (with periodic boundary conditions) are 
considered in order to keep the CPU time into reason- 
able limits. Two types of finite-size effects are present: 
first, the discreteness of the spectrum introduces arti- 
ficial features at low temperatures T < Q.lt and weak 
couphngs [/ < 2i (signaled by (|y)n > 0); second, the 
superconducting phase transition is rounded and corre- 
sponds to the point where the correlation length ^(T) 
becomes larger than the linear system size L. 
Determination ofTc. — Extracted by finite-size analysis 
of very good quality data, the value of Tc is in princi- 
ple free of systematic errors, except a small uncertainty 
(< 5 percents) due to the statistical error and to the fi- 
nite imaginary time discretization At > 0. Given U and 
n, the pair-pair correlation function Ca (Eq.||) is evalu- 
ated for various temperatures T and sizes N. This shows 
clearly that Ca is characterized by a low- and a high- 
temperature regime, related by a transition region that 
becomes sharper and sharper as N increases. The lat- 
ter observation agrees well with the behavior expected in 
the thermodynamic limit, where Ca displays a discon- 
tinuous derivative at the phase transition and becomes 
non-zero only below Tc- This behavior, typical for all the 
parameter values used in our calculations, is shown in 
Fig. H for the special case U = 6t and n = 0.5. Although 
it does not correspond to a genuine phase transition, 
it allows to define a size-dependent transition temper- 
ature Tc{N) which we can use to deduce the value of 
Tc = Tc{oo). A convenient choice for Tc{N) is given by 
the inflexion point of the curve Ca {T, N) versus T ob- 
tained by a (stable) Lagrange polynomial interpolation 
of the DQMC data. Plotting the obtained Tc{N) versus 
N~^, we extrapolate to iV cx) using a linear fit of the 
data, as shown in the inset of Fig. ^ The validity of this 
procedure is confirmed by the evaluation of the specific 
heat cy(r) whose well-defined peak can be used to de- 
fine another size-dependent critical temperature Tl^{N). 
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FIG. 1: Mavn: temperature and size dependences of the pair- 
pair correlation function ^ for the case U = 6t and n = 0.5. 
Inset: extrapolation to the thermodynamic limit of the size- 
dependent critical temperatures Tc{N) and Tl^{N), same U 
and n. 



cv{T) is obtained by the numerical derivative of the ex- 
pectation value of the energy ||l^ . A fit of the values for 
T^*''(7V) with a functional form Tc{oo) = T*^iN)+0{^) 
(corresponding to a superconducting phase transition in 
the universality class of the 3D XY model |l^) is 
shown on the inset of Fig. |l[ It reveals that the finite- 
size corrections to Tc are very weak and in particular 
not larger than the statistical errors resulting from the 
DQMC method. Thus both approaches presented above 
arc fully compatible and yield an uncertainty on the ex- 
trapolated value of Tc which is typically of the order of 
5 percent. 

The critical temperature Tc{U,n). — The above proce- 
dure, applied to a range of parameters U and n, deter- 
mines quantitatively the U-n-Tc phase diagram of the 
AHM. First we consider half-filling (i.e. n = 1) which 
provides a useful check for our method. This case is 
equivalent to the repulsive model that has been recently 
studied by Staudt et al. using the same method § . The 
agreement on the values of TciJJ, n = 1) is almost perfect 
|]lT| ; a small difference (< 3 percents) appearing system- 
atically is due to the extrapolation At — > performed 
by these authors and not done here due to calculation 
time restrictions. Turning now to quarter-filling, we ob- 
tain the results presented on Fig. ||. Before discussing 
the intermediate U regime, we observe that, as long as 
the DMQC method works properly {2t < U < 12t), 
the extreme values of Tc(U) are joining progressively the 
corresponding asymptotic curves, given by the BCS gap 
equation for small U and by the 3D BEC formula for 
large U . Their respective dependences in U follow essen- 
tially Tc (X exp(l//7) and Tc cx 1/C/, with the assumption 
that in the latter case the bosons are noninteracting and 
have an effective hopping amplitude ts = It} jV . In the 
crossover region we observe, as expected, a smooth in- 
terpolation between the BCS and BEC regimes, with a 
maximal value of Tc = 0.35t situated at J7 ~ 8t. It is 
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now interesting to compare our results with those pro- 
posed in recent works. In Fig. ^ the data obtained us- 
ing the Dynamical Mean-Field Theory (DMFT) and a k- 
independent T-matrix approximation (TMA) are also 
plotted, rescaled by a factor ~ 2 so that the dimension- 
less product U times the density of states at the Fermi 
level is the same as in our model ^ . Similarly to the 
half-filling case we observe a good overall agreement 
with DMFT results, the discrepancy at strong coupling 
{U > 6t) being attributed to the mean-field character 
of DMFT; on the other hand TMA clearly fails outside 
the BCS regime. We also mention a recent k-dependent 
T-matrix calculation for U — 4t with a in quanti- 
tative agreement with our results. 

In addition to U, also depends upon the density n. 
Our results show that the function Tc{U = const., n) is 
not monotonic in < n < 1 (1 < n < 2) unlike 
it was previously assumed Q. The maximal transition 
temperature for a given U is situated around n = 0.75 
(1.25). This feature is illustrated on Fig. || for the case 
U = 6t and is reminiscent of the two-dimensional case 
where the higher symmetry of the Hamiltonian at 
half-fiUing (S0(3) instead of U(l)) reduces Tc to zero, as 
discussed recently [0. The appearance of an additional 
charge-density wave (CDW) ordering has been studied 
by means of the corresponding static correlation func- 
tion 

The pairing temperature scale T*{U,n). — As mentioned 
previously, T* is besides Tc another temperature scale 
that characterizes the BCS-BEC crossover. In the case of 
the AHM, T* can be interpreted within a pairing scenario 
as signaling a re-arrangement of fermionic quasiparticles 
into s-wave singlet pairs. As a consequence, the spec- 
tral weight of low-energy spin excitations is reduced and 
the spin response weakens. This process can be studied 
by considering the Pauli susceptibility xp (Eq. Al- 
though T* may not always correspond to a single point, 
but to an extended energy scale, it can nevertheless be 



FIG. 3: Pauli susceptibility xp- Top: T dependence for vari- 
ous values of U (size A'' = 6'^). Bottom: T and A*' dependence 
close to the transition temperature and separation of Tc and 
T* , same symbols as in Fig. |^ (n = 0.5 for both cases). 



identified with the position of the maximum of xp{T) 
JlSf . This definition has the advantage of satisfying the 
expected asymptotic behavior of T*, i.e. T* = Tc in the 
BCS case and T* cx U/\n[U/eF f/^ for the BEC limit 
Q . The way the Pauli spin susceptibility evolves between 
these two regimes is shown on Fig. |[ It is instructive to 
analyze these DQMC results by considering the sensitiv- 
ity of T* to finite-size effects. For the "weak coupling" 
case U < At, one observes that the shape of Xp{T) in 
the region around its maximal value depends strongly on 
the system size N, becoming sharper as N increases. In 
this case the extracted value of T* turns out to be nearly 
equal to Tc, given the accuracy on the numerical results 
(< 5 percents). On the other hand, a "strong couphng" 
behavior appears for U > 5t, characterized by a much 
smoother susceptibility around its maximum. In this re- 
gion finite-size effects have disappeared, indicative of an 
effect characterized by a short coherence length. Here, 
T* is definitely different from Tc. In the interval [Tc, T*] 
the interesting phenomenon of precursor pairing takes 
place, a point which will be further discussed below. We 
can thus present the complete phase diagram on Fig.^ 
by adding the function T*{U,n = const.) that clearly 
displays the two different regimes described above. In 
the weak coupling regime one observes that T* does not 
correspond to a BCS critical temperature extrapolated 
at J7 > 2t. On the strong coupling side T* defines an 
energy scale, which is approximately quantified by the 
errobars on Fig.^ and resembles to a straight line situ- 
ated below the diagonal, in qualitative agreement with 
the asymptotic expression given above. 
Discussion. — A first remark concerns the recent ob- 
servation of a (first-order) phase transition in the non- 
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FIG. 4: U-T phase diagram of the 3D attractive Hubbard 
model at quarter-filhng. The errorbars correspond to the tem- 
perature interval around T* where Xp{T) is less than one per- 
cent smaller than its maximal value, giving thereby an idea 
of the temperature range associated with T* . 

superconducting solution of the AHM 0, ^ . Since this 
state is metastable belovif Tc, it cannot be accessed by 
DQMC (applying a magnetic field would cause a minus- 
sign problem). However the manifestations of this tran- 
sition may be present above as well and the previ- 
ous analysis of the Pauli spin susceptibility constitutes 
an ideal illustration. Indeed, it turns out that the high- 
temperature behavior of xp{T), observed for U < 4t and 
characterized by a monotonic decrease with T, may cor- 
respond to a Fermi liquid normal state where the inter- 
action amounts only to a renormalization of parameters. 
On the other hand, the regime U > 5t, which displays 
the phenomenon of precursor pairing for Tc < T ^ T* , 
fits well to a phase containing "incoherent pairs" |0, |j. 
Consequently a "critical" coupling strength Uc may be 
situated around U = (4.3±0.1)t, as it can be deduced by 
extrapolation at T = on Fig. ^ This value argrees very 
well with the (rescaled) DMFT result 0.56 x W^x 2 « 4.5i, 
W — At being the bandwidth One also remarks 

that Uc does not correspond to the point where the chem- 
ical potential fx (including the Hartree shift —U/2) be- 
comes lower than the bottom of the non-interacting band 
(for n = 0.5, we would get Uc lOt). In fact, to our 
knowledge, there exists no criterion that yields a good 
estimate of Uc in three dimensions. 

In contrast to 3D where the effects of the thermal pair- 
ing fluctuations are rather weak |l^ , in 2D they are 
very important Jl3| leading apparently to a T* joining 
smoothly Tc p9[ . This confirms the observation by Mouk- 
ouri et al. [E0| that precursor phenomena in the AHM 
have two origins: enhanced thermal pairing fluctuations 
(in 2D only) and a strong pairing interaction (in both 
cases). The fact that the AHM contains a transition be- 
tween a Fermi liquid and a state of "incoherent pairs" 
may be of interest in the context of the high- Tc super- 
conductors phase diagram, where the scenario of a hid- 



den quantum phase transition has been proposed pi| . 
Of course the driving parameter in this case is the dop- 
ing and the symmetry of the superconducting phase is 
d-wave. 

The numerical calculations presented in this work were 
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ported by the Swiss National Science Foundation, IR- 
RMA project, the University of Fribourg and the Uni- 
versity of Neuchatel. We thank F. Gebhard, F. Assaad, 
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